Electrical resistivity tomography combined with seismic data estimates heterogeneous distribution of near-seafloor concentrated gas hydrates within gas chimneys

Near-seafloor concentrated gas hydrates (GHs) containing large amounts of methane have been identified at various gas chimney sites. Although understanding the spatial distribution of GHs is fundamental for assessing their dissociation impact on aggravating global warming and resource potential, the spatial distribution of GHs within gas chimneys remains unclear. Here, we estimate the subseafloor distribution of GHs at a gas chimney site in the Japan Sea using marine electrical resistivity tomography data. The resulting two-dimensional subseafloor resistivity structure shows high anomalies (10–100 Ωm) within seismically inferred gas chimneys. As the resistivity anomalies are aligned with high amplitude seismic reflections and core positions recovering GHs, we interpret the resistivity anomalies are near-seafloor concentrated GH deposits. We also detect various distribution patterns of the high resistivity anomalies including 100-m wide and 40-m thick anomaly near the seafloor and 500-m wide anomaly buried 50 m below the seafloor, suggesting that GHs are heterogeneously distributed. Therefore, considering such heterogeneous GH distribution within gas chimneys is critical for in-depth assessments of GH environmental impacts and energy resources.

The resistivity structure obtained from the inversion of the ERT data specified 10-100 Ωm high-resistivity anomalies within seismically inferred gas chimneys.As the resistivity anomalies were aligned with high amplitude seismic reflections and core positions recovering GHs, we interpreted that the resistivity anomalies are nearseafloor concentrated GH deposits.The high resistivity anomalies showed various distribution patterns within the gas chimneys including 100-m wide and 40-m thick anomaly near the seafloor and 500-m wide anomaly buried 50 m below the seafloor.This indicates that GHs and focused gas are heterogeneously distributed within gas chimneys.

Results
The ERT data were obtained from Umitaka Spur along a survey line of approximately 2.5 km 38 (see Method section; Fig. 1b). Figure 2b shows a pseudo-section of the apparent resistivity of observed ERT data, reflecting the seawater and subseafloor resistivity distribution below the survey line.The apparent resistivity also depends on the towing height.Given that the non-GH background resistivity value is 1 Ωm 38 , differences between the observed apparent resistivity data and the 1 Ωm model response provide information on subseafloor high resistivity zones mitigating the effects of the towing height (Fig. 2d,e).Positive differences > 0 Ωm, indicating the existence of high-resistivity anomalies, were observed at a horizontal distance (200-600 m, 900-1200 m and 1200-2300 m).Although a low apparent resistivity zone was observed at a horizontal distance of 900-1200 m, positive differences were observed for source electrodes of C5-COM, C6-COM, and C7-COM in this zone.Therefore, this zone was strongly influenced by seawater owing to the large distance (> 30 m) between the towing system and the seafloor; the low apparent resistivity values do not necessarily indicate the existence of lowresistivity anomalies below the seafloor.Positive differences were observed in seismic blanking zones (Fig. 3a),  38 conducted prior to the ERT survey.The discolored areas comprise microbial mats, carbonate, and gas hydrate (GH).Red and white circles denote the location of piston coring with and without GH samples, respectively 38 .White lines show outlines of GH-bearing areas inferred from seismic surveys; the areas were determined by amplitude anomalies in three-dimensional seismic reflection data and seafloor topography 53 .We generated the maps in (a,b) using the Generic Mapping Tools 5.4.5 (https:// www.gener ic-mappi ng-tools.org/).
indicating that high resistivity zones may exist in the gas chimney structures.Inversion of the observed ERT data facilitates the incorporation of seawater and towing height effects and quantitatively determines the resistivity structure below the seafloor.
The inversion model (Fig. 3b) revealed the prominent high-resistivity zones of R1-R6 (> 10 Ωm) with various distribution patterns.The recovery of the high-resistivity zones R1-R6, can be explained by positive differences between the observed data and the 1 Ωm seafloor model response (Fig. 2e).Seismic blanking zones were   observed at and below R1-R6, indicating R1-R6 are related to the gas chimney structures (Fig. 3a), whereas non-GH background resistivity zones (~ 1 Ωm) was in the non-gas chimney regions (regions without seismic blanking).R1 to R5 were located near the seafloor, whereas R6 was buried at 50 mbsf.The near-seafloor resistors R1-R5 differed in their horizontal reach and vertical thickness.For example, R1 was 70 m wide and 100 m thick, whereas R2 was 100 m wide and 40 m thick.R1 and R4 extended to the deeper subseafloor regions, whereas R2, R3, and R5 extended to tens of meters below the seafloor.R2 and R5 coincided with the positions of the piston cores, which recovered GHs 38 .R1, R2, and R5 were on the mound structure, whereas R3 and R4 were on the  38 .Green arrows show GH areas inferred from seismic surveys 53 (Fig. 1b).Black arrows show discolored seafloor areas detected based on video observations 38 .Orange arrows show methane plumes zones 21 .
pockmark structure; R1-R5 were consistent with GH-bearing areas inferred from seismic surveys 53 .Methane plumes were observed at R2 and R5 21 .R2-R5 were located at zones of near-seafloor high amplitude reflections in the seismic section (Fig. 3a,b).Note that the three-dimensional (3D) seismic survey was aimed at deep reservoirs at 4000-5000 mbsf and was not intended for shallow GH; the wavelength was ~ 50 m for 30-40 Hz center frequency and 1800 m/s velocity 59 .
Subseafloor GH/gas saturation in pore spaces (S g ) was calculated from the resistivity model using Archie's equation: a = 1.0, n = 2.0, and m = 2.5 (see Method section; Fig. 3c); the porosity values are described in Eq. ( 6).High saturation rates > 60% were estimated in the high-resistivity zones of R1-R6 (Fig. 3c).The high GH saturation rates for R1-R6 were also computed with different m = 2.2-2.8.S g of 100 Ωm (R2) calculated with m = 2.2, 2.5, and 2.8 were 91.4%, 91.0%, and 90.5%, respectively.S g of 13 Ωm (R5) calculated with m = 2.2, 2.5, and 2.8 were 76.2%, 75.0%, and 73.7%, respectively.These results indicated that the S g estimation of R1-R6 only differed by a maximum of 3% with an m = 2.2-2.8 and the influence of the cementation factor on the saturation estimate decreased when the formation resistivity or saturation was higher.

Discussion
Our resistivity model reveals near-seafloor high-resistivity zones of R1-R5, interpreted as near-seafloor concentrated GHs.The obtained resistivity values for R1-R5 (10-100 Ωm) are consistent with those reported for GHs at other gas chimney sites 34,45 .R1-R5 are located on seafloor mounds and pockmarks; the distribution of R1-R5 is consistent with GH-bearing areas inferred from seismic surveys 53 .R1, R2, R4, and R5 overlap with discolored areas detected based on seafloor video observations 38 ; the discolored areas consist of microbial mats, carbonate, and GH.R2 and R5 coincide with the positions where GHs were identified by piston cores.The piston cores collected massive-type GH samples at R2 38 .The piston core at R5 recovered the soupy mud, implying the dissociation of GH during piston core recovery 38 .The higher resistivity values of R2 than R5 are consistent with the core recovery of GHs.High amplitude reflections observed at R2-R5, suggest that R2-R5 are near-seafloor concentrated GHs.High amplitude reflections are not observed at R1, implying that GH saturation rates of R1 are lower than those of R2-R5.Seismic blanking zones below R1-R5 imply high-flux gas upwells to R1-R5.Although resistivity alone cannot discriminate GH or gas 34 , our resistivity model combined with preceding seismic data and core recovery observations reinforces our interpretation that GH is the predominant constituent for the high-resistivity zones of R2-R5, with high-flux gas inputs resulting in GH accumulation at R1-R5 (Fig. 4).
High GH/gas saturations exceeding 60% in R1-R4 indicate the presence of high saturation GH deposits with substantial volumes.These deposits exhibit horizontal dimensions ranging from 50 to 150 m and vertical dimensions spanning 40 m to 80 m.In the Joetsu Basin, which encompasses our study area, coring was performed at nine drill sites within seismically inferred gas chimneys and core analysis reported that the average volume fraction of GHs at each drill site ranges from 35% to 86% of the sedimentary sequences 20 .Comparing these values with our GH saturation model suggests that our GH saturation model provides reasonable estimates.
Horizontally elongated R6, buried at a depth of 50 mbsf, differs geometrically from the near-seafloor GHs of R1-R5.Such horizontally elongated GHs have been found above BSRs in various GH areas.In facts, BSRs are observed below R6 at a horizontal distance of 1300-1500 m (Fig. 3a).Our ERT data cannot constrain the bottom of R6 due to the shallow exploration depth; however, the bottom of the GH-bearing zone of R6 is constrained by the BGHS depth or BSRs.The BGHS depth in this area is reportedly 115 mbsf 21 .The thickness of the GH-bearing zone of R6 will be 70 m if it extends to the BGHS depth, corresponding to the upper-bound estimate of the GHbearing zone of R6.Moreover, given that high amplitude reflections are not observed at R6, GH saturation rates of R6 are lower than those of R2-R5, while the free gas saturation rates of R6 are higher than those of R2-R5.
The GH/gas saturation rates of R2 and R3 are the highest near the seafloor and decrease with depth.A thermodynamic equilibrium model, which considers a three-phase zone, can explain much higher GH saturation near the seafloor than in the deep parts of a gas chimney 4,60,61 .That is, free gas upwelling from depths enters the GH stability zone, where the reaction of free gas with seawater forms GH; porewater salinity increases as a result of GH formation; increasing salinity brings the local thermodynamic condition closer to a three-phase equilibrium of seawater, GH, and free gas.GH formation continues until the system reaches a three-phase equilibrium due to the sufficiently elevated salinity.The free gas can rise to the seafloor after the GH stability zone reaches the three-phase equilibrium.In the case of high gas flux with no water flux in sand, high salinity, low water content, and high GH saturation are required to maintain the three-phase equilibrium conditions near the seafloor 61 .In contrast, low salinity, high water content, and low GH saturation explain the three-phase equilibrium conditions in deeper parts, including immediately above the BGHS 61 .These numerical predictions are consistent with the GH saturations of R2 and R3 observed in this study.The conductor CD1 below R2 exhibits unusually high conductivities within the GH stability zone, suggesting that it might be an artifact resulting from the inversion.As CD1 is situated in a circular GH area inferred from seismic surveys 53 (Figs.1b and 3b), 3D effects of the circular structure to our 2D inversion could lead to the artifact.
Our resistivity model estimates potential GH zones within the seismically inferred gas chimneys by combining the seismic data, the piston cores, and the numerical thermodynamic model.The estimated GH saturation rate for R2 and R3 is the highest near the seafloor and decreases with depth; this is explained well by the thermodynamic equilibrium model 60 .Various distribution patterns of potential GH deposits are found within the seismically inferred gas chimneys: (1) GHs extending from the seafloor to deeper zones (R1 and R4), (2) GHs extending from the seafloor to several tens of meters below the seafloor (R2, R3, and R5), and (3) horizontally elongated GH buried at a depth of 50 mbsf (R6).This indicates that GHs and focused gas flows are heterogeneously distributed within gas chimneys.Such a variety of GH distributions were not considered when determining their dissociation impact on environments.Incorporating the near-seafloor GH distribution into analyses of their dissociation impact on environments will improve the assessment accuracy.Moreover, considering the various GH distributions will improve the accuracy of GH resource estimations.

Deep-towed marine ERT method
Deep-towed marine ERT methods used electric current sources and receivers to investigate subseafloor resistivity structures 38,48,62 .The ERT data used in this study were acquired by Goto et al. 38 during the research cruise KY05-08 in August 2005.The ERT system comprised eight source and two receiver electrodes attached to a 160-m long floating cable (Fig. 2a) that was behind the main system, including cameras, conductivity-temperature-depth (CTD) sensors, an altimeter, and an alternating current transformer 38 .The ERT system was towed by a ship at an altitude of less than a few tens of meters from the seafloor with a towing speed of ~ 1 knot (0.5 m/s).These short-spacing electrode arrays with the near-seafloor-towed ERT system enabled near-seafloor resistivity imaging.The ERT system used a source signal comprising sinusoidal with a period of 4 s; a source signal was sequentially sent to the seven dipoles for one period at a constant voltage of 72 V peak to peak.Accordingly, one source signal sequence took 28 s.The ERT system recorded the electrical potential at the P1 and P2 electrodes.The sampling rate of source current and receiver potential data was 0.2 s.Apparent resistivity was computed by stacking data for three sequences, i.e., 84 s, corresponding to ~ 50 m array movement.To mitigate the effects of the metallic cable and deep-tow frame of the ERT system, the apparent resistivity data were calibrated by conducting measurements in the water column at a depth of 400 m.The calibration factors of the apparent resistivity were 1.05, 1.03, 1.02, 1.14, 1.13, 1.10, and 1.08 for source dipoles of C1-COM, C2-COM, C3-COM, C4-COM, C5-COM, C6-COM, and C7-COM, respectively.An acoustic supershort-baseline (SSBL) navigation system provided positioning information on the transponders attached to the deep-tow frame and the cable tail.Goto et al. 38 provide further details on the ERT system.
The apparent resistivity formulas help interpret ERT data features 38,48,62 .The apparent resistivity was determined from the observed voltage, source current, and geometries of the source and receiver electrodes using Eq.(1): where ρ a is apparent resistivity (Ωm), V is receiver voltage (V), I is source current (A), and r is the distance between electrodes (m).Subscripts c and p indicate the current and receiver electrodes, respectively; ρ a represents the resistivity, assuming a homogeneous resistivity structure.

Inversion of observed ERT data
The ERT inversion estimated the subseafloor resistivity structure from the apparent resistivity data.We used the 2D inversion code developed by Ishizu et al. 48to obtain the 2D resistivity structure below the towing profile.The inversion minimizes the following objective function U, where m is a model vector, m 0 is a prior model vector, d is the observed data vector, F[m] represents the forward modeling response, C m represents a model covariance matrix, C d is a data covariance matrix, and µ is a tradeoff parameter between data misfit and model roughness.The first derivative roughness penalty was used for C −1 m .Owing to the non-linearity of the ERT inversion, the inversion applied an iterative approach based on the Gauss-Newton algorithm.Inversion with the Gauss-Newton algorithm generally converges after a small (1) iteration number 63 .The tradeoff parameter µ was determined by an L-curve criterion 64 .The root mean square (RMS) data misfit was defined to measure overall data fitting per Eq. ( 3): where N is the number of data.The forward modeling used the finite element method with unstructured triangular meshes because the unstructured triangular meshes can model the seafloor topography efficiently 48 .The total number of apparent resistivity values input into the inversion code (N) was 357.The initial and prior models for the inversion consisted of 1.0 Ωm homogeneous subseafloor and 0.35 Ωm seawater layers.This sea resistivity value was the average value of the seawater resistivity measured using the CTD sensor during the ERT survey 38 .The mesh number for forward modeling was 20,318.For accurate computation, each mesh size was set to < 30 m 2 in areas with high data sensitivity (i.e., a horizontal distance of −100 to 2600 m and a depth of from 800 to 1040 m).The resistivity of the sea layer was fixed at 0.35 Ωm during the inversion, which only updated the subseafloor region.Consequently, the total number of unknown model parameters was 9935.A minimum error setting of 2% was used for inversion.Seafloor topography was modeled as the sum of the CTD depth and altitude from the seafloor, measured using an altimeter.The altimeter data were missing for the seafloor topography at a horizontal distance of 1000-1250 m; thus, the topography of this section was supplemented with 50 m × 50 m bathymetric data (Fig. 1b).The 2D inversion, with a fixed tradeoff parameter µ = 30, yielded a resistivity model below the seafloor after three iterations (Fig. 3b), and the RMS data misfit was 1.76.A tradeoff parameter of 30 was selected as the curvature of the L-curve was the largest for the parameter (Supplementary Fig. 1).The model explains the major features of the observed data (Fig. 2b,c; Supplementary Fig. 2).Meanwhile, the data fit worsened at 250-300 m and 1200-1250 m due to possible 3D topography effects.

Error of tail angle and towing altitude
The ERT data included the systematic errors of the tail angle and towing altitude, which can distort inversion results.The standard error of the tail depth monitored by the SSBL system was ~ 0.5 m 38 ; thus, the tail angle error was < 0.2°.A synthetic model study demonstrated that a tail error up to 1° did not strongly influence the inversion result of 2D ERT data 65 .Thus, we concluded that the systematic error of the tail angle did not significantly impact the inverted model.
The towing altitudes were measured using ultrasonic reflections from the seafloor 38 and may have been underestimated as complex structures, such as mounds, can cause wave reflections from the side rather than directly below the towing system.Moreover, altitude data were missing at horizontal distance of 1000-1250 m and the altitudes may have been largely overestimated or underestimated.A synthetic model study reported that if the towing height was underestimated by 2.5 m (e.g., recorded height was 7.5 m, the actual towing height was 10 m), the inverted resistivity model did not differ significantly 65 .The error of the recorded altitude data was considered to be < 2.5 m.Although R3 and R4 were located on the altitude missing positions, the seismic section specified near-seafloor high amplitude reflections at R3 and R4 positions; this consistency supports the credibility of the resistivity model.

Estimation of GH/gas saturation
The GH/gas saturation rates were estimated from the subseafloor resistivity values using Archie's law-an empirical law that describes the relationship between electrical resistivity and pore spaces in sedimentary rocks 66 .Saturation estimates obtained from resistivity data using Archie's law have been used in various studies [34][35][36]67 . Arcie's law is expressed by Eq. ( 4): where S w is water saturation rate in pore spaces, R w is the resistivity of pore water (Ωm) set to the seawater resistivity of 0.35 Ωm, R t represents bulk resistivity (Ωm), φ is the porosity of sediment, a is the tortuosity factor, m is the cementation factor, and n is the saturation exponent 36,66 .Assuming that GH and/or gas fills the remaining pore space, the GH/gas saturation rate in the pore spaces (S g ) is defined by Eq. ( 5): The tortuosity factor (a) and saturation exponent (n) were selected as 1.0 and 2.0, respectively, according to a study by Schwalenberg et al. 34 ; a cementation factor (m) ranging from 2.2 to 2.8 has been used to estimate saturation in drilled GH areas 34 .Thus, given that the information required to determine the value of m was unavailable for this study, we applied the m = 2.2-2.8range.Goto et al. 68 reported a vertical porosity function for the Joetsu Basin, including our study area Umitaka Spur, using Eq. ( 6): where z is the depth below the seafloor (m); this function was derived from core samples obtained from a depth up to 3081 mbsf.The closest drilling site (MD179-3304) was within 1 km of the ERT profile, and the geological setting of the drill sites was comparable to our survey area.We applied Eq. ( 6) to the porosity values used in Archie's law to obtain S w and S g .

Figure 1 .
Figure 1.Map of the study area.(a) Location of Umitaka Spur, the eastern margin of the Japan Sea, is indicated by a yellow star.(b) Event map of Umitaka Spur.The magenta line represents a survey profile with a towed electrical resistivity tomography (ERT) system 38 .Black lines indicate discolored seafloor areas, detected based on video observations along a towing profile38 conducted prior to the ERT survey.The discolored areas comprise microbial mats, carbonate, and gas hydrate (GH).Red and white circles denote the location of piston coring with and without GH samples, respectively38 .White lines show outlines of GH-bearing areas inferred from seismic surveys; the areas were determined by amplitude anomalies in three-dimensional seismic reflection data and seafloor topography53 .We generated the maps in (a,b) using the Generic Mapping Tools 5.4.5 (https:// www.gener ic-mappi ng-tools.org/).
Difference between the observed data and response from the initial

Figure 2 .
Figure 2. Electrode configuration of the electrical resistivity tomography (ERT) system, observed ERT data, and model responses.(a) Electrode location of the ERT system 38 .The ERT system uses eight current electrodes, C1-C7 and COM, and two receiver electrodes, P1 and P2.Current was transmitted using electrode pairs: C1-COM, C2-COM, C3-COM, C4-COM, C5-COM, C6-COM, and C7-COM.(b) Observed apparent resistivity data 38 .(c) Response from the inverted resistivity model shown in Fig. 3b.(d) Response from an initial model comprising a 1.0 Ωm homogeneous subseafloor layer and a 0.35 Ωm seawater layer.(e) Difference between the observed data and the initial model response.(d) and (e) highlight data that the homogeneous sub-seafloor structure of 1 Ωm cannot explain.The vertical axis is based on source electrode pairs.The C7-COM pair reflects deeper resistivity structures than C1-COM because of longer source-receiver septation.Data are shown with southwest on the left (start point of towing) and northeast on the right (end point of towing).The observed apparent resistivity at a horizontal distance of 900-1200 m is low because of the high tow height of the system.

Figure 3 .
Figure 3. Seismic section, inverted resistivity model, and subseafloor saturation rates.(a) Seismic section below the electrical resistivity tomography (ERT) towing profile (Fig. 1b).The seismic section is extracted from threedimensional seismic survey data 59 .Two-way travel time of 0.1 s corresponds to a depth of 85-90 m for velocities of 1700-1800 m/s 59 .A magenta rectangular represents an outline of the resistivity model in (b).Green and light blue dashed circles show near-seafloor high amplitude reflections and bottom simulating reflectors (BSRs), respectively.Near-vertical, dark-blue dashed curves indicate the boundaries of seismic blanking zones.(b) Vertical cross-section of the inverted resistivity model below the towing profile.R1-R6 and CD1 are high (> 10 Ωm) and low (< 1.0 Ωm) resistivity zones, respectively.(c) Distribution of gas hydrate (GH)/gas saturation rates.Saturation rates are derived from the subseafloor resistivity values in (b) using Archie's law.Small black squares and triangles indicate the locations of the head and tail of the towing system of ERT data.Each location of the head and tail is connected by lines.Red and white circles indicate the positions of piston cores with and without GH, respectively, within 150 m of the ERT towing profile38 .Green arrows show GH areas inferred from seismic surveys53 (Fig.1b).Black arrows show discolored seafloor areas detected based on video observations38 .Orange arrows show methane plumes zones21 .

Figure 4 .
Figure 4. Schematic diagram of the gas hydrate (GH) evolution model at a gas chimney site in the Japan Sea.The white zones are areas with peak GH/gas saturation > 60% and red arrows represent focused gas pathways.Light green lines indicate uncertain bottoms of the GHs that are unconstrained by our electrical resistivity tomography data.Blue lines and yellow dashed curves represent bottom simulating reflectors (BSRs) and outlines of seismically inferred gas chimneys.The figure is not to scale.